Подгружаем датасет
insurance <- read_csv("./insurance_cost.csv")
## Rows: 1338 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): sex, smoker, region
## dbl (4): age, bmi, children, charges
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Ради большего интереса переименуем его.
endurance <- insurance
What endurance looks like
Смотрим на переменные
str(endurance)
## spc_tbl_ [1,338 × 7] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ age : num [1:1338] 19 18 28 33 32 31 46 37 37 60 ...
## $ sex : chr [1:1338] "female" "male" "male" "male" ...
## $ bmi : num [1:1338] 27.9 33.8 33 22.7 28.9 ...
## $ children: num [1:1338] 0 1 3 0 0 0 1 3 2 0 ...
## $ smoker : chr [1:1338] "yes" "no" "no" "no" ...
## $ region : chr [1:1338] "southwest" "southeast" "southeast" "northwest" ...
## $ charges : num [1:1338] 16885 1726 4449 21984 3867 ...
## - attr(*, "spec")=
## .. cols(
## .. age = col_double(),
## .. sex = col_character(),
## .. bmi = col_double(),
## .. children = col_double(),
## .. smoker = col_character(),
## .. region = col_character(),
## .. charges = col_double()
## .. )
## - attr(*, "problems")=<externalptr>
summary(endurance)
## age sex bmi children
## Min. :18.00 Length:1338 Min. :15.96 Min. :0.000
## 1st Qu.:27.00 Class :character 1st Qu.:26.30 1st Qu.:0.000
## Median :39.00 Mode :character Median :30.40 Median :1.000
## Mean :39.21 Mean :30.66 Mean :1.095
## 3rd Qu.:51.00 3rd Qu.:34.69 3rd Qu.:2.000
## Max. :64.00 Max. :53.13 Max. :5.000
## smoker region charges
## Length:1338 Length:1338 Min. : 1122
## Class :character Class :character 1st Qu.: 4740
## Mode :character Mode :character Median : 9382
## Mean :13270
## 3rd Qu.:16640
## Max. :63770
Все выглядит пристойно.
Оценим пропущенные значения.
endurance %>%
is.na() %>%
sum()
## [1] 0
Их нет, и это прекрасно.
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
В данных не было указано, в какой валюте осуществлялись страховые выплаты. Выберем в качестве валюты наиболее вероятный вариант – монгольские тугрики.
charges_density <- ggplot(endurance)+
geom_density(aes(x=charges),
fill="white",
color="green",
alpha=0.9)+
theme_dark()+ # Помню, что обычно рекомендуют theme_bw(), но раз в задании рекомендация поиграть с темами -- пусть будет соответствовать названию датасета
geom_vline(aes(xintercept=mean(charges)),
color="red",
size=1.5,
linetype="dashed")+
geom_vline(aes(xintercept=median(charges)),
color="blue",
size=1.5,
linetype="dashed")+
annotate(geom = "text",
x = 40000,
y = 0.00005,
label = glue("Среднее для страховых выплат составляет\n{round(mean(endurance$charges),0)} монгольских тугриков"),
color = "red",
size=5)+
annotate(geom = "text",
x = 40000,
y = 0.00003,
label = glue("Медиана для страховых выплат составляет\n{round(median(endurance$charges),0)} монгольских тугриков"),
color = "blue",
size=5)+
labs(x = "Страховые выплаты",
y = "Плотность вероятности",
title = "Рис. 2 График плотности вероятности\nстраховых выплат")+
theme(axis.text.y = element_text(size=16),
axis.title.x = element_text(size=19),
axis.title.y = element_text(size=19),
plot.title = element_text(size=20, hjust=0.5),
legend.title = element_text(size=21),
legend.text = element_text(size=18))
charges_density
# Функция для массового создания боксплотов в чанке setup
endurance_boxps <- lapply(endurance %>% select(where(is.character)) %>% colnames, function(x) boxpCreator(endurance, x, "charges"))
endurance_boxps
## [[1]]
##
## [[2]]
##
## [[3]]
boxplots <- ggarrange(endurance_boxps[[1]],
endurance_boxps[[2]],
endurance_boxps[[3]],
ncol=3, nrow=1)
combo_plot <- ggarrange(charges_density, boxplots,
ncol=1, nrow=2)
annotate_figure(combo_plot, top = text_grob("Рис.4 Вероятность страховых выплат различного размера\nи распределение страховых выплат в зависимости от пола, курения и региона проживания\n",
color = "black", face = "italic", size = 25))
Базовый вариант
charges_density+
facet_grid("region")+
labs(title = "Рис. 5 График плотности страховых выплат по регионам")
age_charges_scatter <- ggplot(endurance)+
geom_point(aes(age, charges), colour="violet")+
theme_bw()+
labs(x = "Возраст",
y = "Сумма страховки",
title = "Рис. 6 Диаграмма рассеяния\nвозраст-страховые выплаты")+
theme(axis.text.x = element_text(size=14), # Как заказывали
axis.text.y = element_text(size=10),
axis.title.x = element_text(size=10),
axis.title.y = element_text(size=10),
plot.title = element_text(size=20, hjust=0.5),
legend.title = element_text(size=21),
legend.text = element_text(size=18))
age_charges_scatter
age_charges_scatter+
geom_smooth(aes(x=age, y=charges),
color = "green",
linewidth = 1.5,
method = lm,
se = FALSE)
## `geom_smooth()` using formula = 'y ~ x'
age_charges_scatter+
geom_point(aes(age, charges, colour=smoker))+
geom_smooth(aes(age, charges, linewidth=smoker, group=smoker),
method=lm,
se=FALSE)+
labs(colour="Курение",
size="Курение")
## `geom_smooth()` using formula = 'y ~ x'
ggplot(endurance, aes(bmi, charges, colour=smoker))+
geom_point()+
geom_smooth(aes(group=smoker),
method=lm,
se=FALSE)+
theme_bw()+
labs(x = "bmi",
y = "Сумма страховки",
title = "Рис. 7 Диаграмма рассеяния\nИМТ-страховые выплаты",
colour="Курение")+
theme(axis.text.x = element_text(size=18),
axis.text.y = element_text(size=18),
axis.title.x = element_text(size=20),
axis.title.y = element_text(size=20),
plot.title = element_text(size=20, hjust=0.5),
legend.title = element_text(size=21),
legend.text = element_text(size=18))
Я бы сказал, что график по форме напоминает два шашлыка на шампуре под углом 35 градусов друг к другу. Из литературных данных известно, что подобная shashlyk-shaped форма графика наглядно показывает, что курение, вероятно, намного сильнее взаимосвязано с суммой страховых расходов, чем индекс массы тела.
Вопрос не про страховку: отличается ли ИМТ у имеющих одного ребенка некурящих мужчин и женщин.
question1 <- endurance %>%
filter(children == 1 & smoker == "no") %>%
ggplot()+
geom_boxplot(aes(x=sex, y=bmi, fill=sex))+
theme_bw()+
labs(x = "Пол",
y = "ИМТ",
title = "Рис. 7 Распределение ИМТ по полу у некурящих\nлиц, имеющих одного ребенка",
fill = "Пол")+
theme(axis.text.x = element_text(size=18),
axis.text.y = element_text(size=18),
axis.title.x = element_text(size=20),
axis.title.y = element_text(size=20),
plot.title = element_text(size=20, hjust=0.5),
legend.title = element_text(size=21),
legend.text = element_text(size=18))
question1
Визуально кажется, что нет. Проверим.
endurance %>%
filter(children == 1 & smoker == "no") %>%
group_by(sex) %>%
do(broom::tidy(shapiro.test(.$bmi)))
## # A tibble: 2 × 4
## # Groups: sex [2]
## sex statistic p.value method
## <chr> <dbl> <dbl> <chr>
## 1 female 0.994 0.837 Shapiro-Wilk normality test
## 2 male 0.981 0.0661 Shapiro-Wilk normality test
Не можем отвергнуть нулевую гипотезу о том, что распределения ИМТ в целевых группах нормальны. Используем тест Стьюдента для проверки гипотезы о равенстве средних.
stat.test <- endurance %>%
filter(children == 1 & smoker == "no") %>%
t.test(bmi ~ sex, data = .) #%>%
#add_xy_position(x = "sex")
stat.test
##
## Welch Two Sample t-test
##
## data: bmi by sex
## t = -0.95426, df = 256.89, p-value = 0.3408
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
## -2.2035331 0.7650241
## sample estimates:
## mean in group female mean in group male
## 30.20936 30.92862
Нулевую гипотезу отвергнуть не можем. Средние по нашим данным не различаются.
Имеем одну категориальную и одну непрерывную числовую переменную.
Оптимальным способом отображения данных такого типа являются графики
типов боксплот, скрипкаграфик вайолинплот и джиттерплот
Как связаны между собой возраст, размер страховых выплат и ИМТ у курящих и некурящих мужчин?
plot_ly(
endurance %>% filter(sex == "male"),
x = ~ age,
y = ~ bmi,
z = ~ charges,
size = 1,
color = ~ smoker,
colors = c("#FF0033", "#00CC00"),
alpha = 0.5,
width = 1.5,
height = 1.5
) %>%
layout(
title = "Распределение возраста, ИМТи размера выплаченной страховки\nу курящих и некурящих мужчин",
xaxis = list(title = "Возраст",
zeroline = FALSE),
yaxis = list(title = "Индекс массы тела",
zeroline = FALSE)
)
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
Вырисовываются три паралельные группы, в которых с увеличением возраста растет и сумма выплаченной страховки (как у курящих, так и у некурящих, причем у последних суммы были заметно ниже). В случае с индексом массы тела его увеличение приводило к увеличению страховки только у некурящих пациентов.
Имеем три типа переменных. Посмотреть на их взаимное распределение удобно на диаграмме рассеяния (скаттерплоте). Ну а визуализация в плотли выбрана потому, что это забавно выглядит (этот пункт был выполнен до публикации второй домашки, а после передлывать как-то не пришлось по душе).
Как предсказать размер страховки, исходя из имеющихся данных?
Попробуем поиграть с линейной регрессией.
# Стандартизуем количественные переменные
endurance_scaled <- endurance %>%
mutate(across(where(is.numeric), scale))
# Построим модель
charge_lm <- lm(charges ~ bmi + age + children + sex + smoker + region, data = endurance_scaled)
summary(charge_lm)
##
## Call:
## lm(formula = charges ~ bmi + age + children + sex + smoker +
## region, data = endurance_scaled)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.9335 -0.2352 -0.0811 0.1151 2.4767
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.34822 0.03189 -10.920 < 0.0000000000000002 ***
## bmi 0.17081 0.01440 11.860 < 0.0000000000000002 ***
## age 0.29800 0.01380 21.587 < 0.0000000000000002 ***
## children 0.04733 0.01372 3.451 0.000577 ***
## sexmale -0.01084 0.02749 -0.394 0.693348
## smokeryes 1.96932 0.03412 57.723 < 0.0000000000000002 ***
## regionnorthwest -0.02915 0.03933 -0.741 0.458769
## regionsoutheast -0.08547 0.03953 -2.162 0.030782 *
## regionsouthwest -0.07928 0.03947 -2.009 0.044765 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5006 on 1329 degrees of freedom
## Multiple R-squared: 0.7509, Adjusted R-squared: 0.7494
## F-statistic: 500.8 on 8 and 1329 DF, p-value: < 0.00000000000000022
Выглядит неплохо уже на текущем этапе –
поправленный R-квадрат = 0.75
Поищем мультиколлинеарность в данных (вот здесь мы и применим визуализацию)
endurance_corr <- corr.test(endurance %>% select(where(is.numeric)& !charges), method="spearman")
corr_plot <- ggcorrplot(endurance_corr$r,
type = "lower",
outline.col = "white",
lab=TRUE,
method = "circle",
p.mat = endurance_corr$p,
title = "Скоррелированность количественных данных\nв датасете по страховке",
legend.title = "Коэффициент\nкорреляции")
corr_plot+
labs()
На первый взгляд коллинеарности нет. Проверим формальными тестами.
vif(charge_lm)
## GVIF Df GVIF^(1/(2*Df))
## bmi 1.106630 1 1.051965
## age 1.016822 1 1.008376
## children 1.004011 1 1.002003
## sex 1.008900 1 1.004440
## smoker 1.012074 1 1.006019
## region 1.098893 3 1.015841
Все действительно в порядке.
Сделаем шаг backward selection:
drop1(charge_lm, test = "F")
## Single term deletions
##
## Model:
## charges ~ bmi + age + children + sex + smoker + region
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 333.03 -1842.76
## bmi 1 35.25 368.28 -1710.15 140.6627 < 0.00000000000000022 ***
## age 1 116.77 449.80 -1442.60 465.9837 < 0.00000000000000022 ***
## children 1 2.98 336.01 -1832.82 11.9063 0.000577 ***
## sex 1 0.04 333.07 -1844.60 0.1556 0.693348
## smoker 1 834.95 1167.98 -165.84 3331.9680 < 0.00000000000000022 ***
## region 3 1.59 334.62 -1842.38 2.1173 0.096221 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Удалим пол, как наименее значимый предиктор.
charge_lm2 <- update(charge_lm, .~. - sex)
drop1(charge_lm2, test = "F")
## Single term deletions
##
## Model:
## charges ~ bmi + age + children + smoker + region
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 333.07 -1844.60
## bmi 1 35.22 368.28 -1712.12 140.6226 < 0.00000000000000022 ***
## age 1 116.95 450.02 -1443.95 466.9969 < 0.00000000000000022 ***
## children 1 2.97 336.04 -1834.71 11.8706 0.000588 ***
## smoker 1 838.82 1171.89 -163.37 3349.5461 < 0.00000000000000022 ***
## region 3 1.59 334.66 -1844.23 2.1166 0.096315 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Еще шаг. Теперь удаляем регион.
charge_lm3 <- update(charge_lm2, .~. - region)
drop1(charge_lm3, test = "F")
## Single term deletions
##
## Model:
## charges ~ bmi + age + children + smoker
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 334.66 -1844.23
## bmi 1 34.70 369.36 -1714.24 138.203 < 0.00000000000000022 ***
## age 1 117.94 452.60 -1442.28 469.789 < 0.00000000000000022 ***
## children 1 2.96 337.62 -1834.43 11.809 0.0006077 ***
## smoker 1 841.77 1176.43 -164.19 3352.911 < 0.00000000000000022 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(charge_lm3)
##
## Call:
## lm(formula = charges ~ bmi + age + children + smoker, data = endurance_scaled)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.98249 -0.24119 -0.08147 0.11497 2.43679
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.40266 0.01536 -26.211 < 0.0000000000000002 ***
## bmi 0.16207 0.01379 11.756 < 0.0000000000000002 ***
## age 0.29916 0.01380 21.675 < 0.0000000000000002 ***
## children 0.04713 0.01372 3.436 0.000608 ***
## smokeryes 1.96626 0.03396 57.904 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5011 on 1333 degrees of freedom
## Multiple R-squared: 0.7497, Adjusted R-squared: 0.7489
## F-statistic: 998.1 on 4 and 1333 DF, p-value: < 0.00000000000000022
Наш R-квадрат почти не изменился, зато теперь нам достаточно четырех
переменных вместо шести (что весьма важно – мы исключили достаточно
трудно интерпретируемую в рамках линейной регрессию переменную
region).
Коррплот сделан как раз для того, чтобы искать скоррелированные переменные. В нашем случае его наглядность может несколько страдать из-за малого количества числовых переменных, одако, когда у нас 10-15 оных, коррплоты могут очень облегчить поиск коллинеарных переменных.
У целевого графика есть расхождения с подписями. Если делать его с
разбиением на группы в соответствии с легендой, то точек на скаттерплоте
меньше, чем на целевом. Если же сбросить нижнюю границу возраста в
первой группе (что приведет к включению лиц в возрасте 18-20 лет в
группу 21-34), график принимает необходимый вид. Поскольку
целью задания является именно построение идентичного графика, я занулил
нижнюю границу в первой возрастной группе.
endurance %>%
mutate(age_group=as.factor(case_when(
age < 35 ~ "age: 21-34",
age >= 35 & age < 50 ~ "age: 35-49",
age >= 50 ~ "age: 50+"
))) %>%
na.omit() %>%
ggplot(aes(x=bmi, y=log(charges)))+
geom_point(color = "#330066", alpha = 0.5, size=2)+
facet_grid(cols=vars(age_group))+
geom_smooth(aes(colour=age_group), method=lm)+
theme_minimal()+
labs(title="Отношение индекса массы тела к логарифму трат по возрастным группам")+
theme(
axis.text.x = element_text(size=16),
axis.text.y = element_text(size=16),
axis.title.x = element_text(size=19),
axis.title.y = element_text(size=19),
plot.title = element_text(size=23, hjust=5),
legend.title = element_text(size=19),
legend.text = element_text(size=16),
legend.position = "bottom",
legend.key.size = unit(1.5, "line"),
strip.text.x = element_text(size=16)
)
## `geom_smooth()` using formula = 'y ~ x'
Смею тешить себя надеждой, что